rm(list =ls())
options(scipen=999)
gc()
packages <-c("tidyverse","estimatr","plm","stargazer",
             "fastDummies","ICCbin","ihs","readstata13","xtable")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)

setwd("PUT YOUR DIRECTORY HERE")

# load data
data <- read.csv("./Datasets/us_generalizability.csv")

# rearrange data
data <- data%>%
  dplyr::mutate(NSZZ_members_1981_capita = NSZZ_members_1981/Pop_1980,
                Investmens_1976_80_capita = Investments_1976_80/((Pop_1975 + Pop_1980)/2),
                Employed_industry_1975_capita = Employed_industry_1975/Pop_1975,
                Employed_industry_1980_capita = Employed_industry_1980/Pop_1980,
                Employed_agricult_1975_capita = Employed_agricult_1975/Pop_1975,
                Employed_agricult_1980_capita = Employed_agricult_1980/Pop_1980,
                Migr_ratio_urban_1975_capita = Migr_ratio_urban_1975/Pop_1975,
                Migr_ratio_urban_1980_capita = Migr_ratio_urban_1980/Pop_1980,
                Migr_ratio_rural_1975_capita = Migr_ratio_rural_1975/Pop_1975,
                Migr_ratio_rural_1980_capita = Migr_ratio_rural_1980/Pop_1980,
                Urban_pop_1975_capita = Urban_pop_1975/Pop_1975,
                Urban_pop_1980_capita = Urban_pop_1980/Pop_1980,
                Rural_pop_1975_capita = Rural_pop_1975/Pop_1975,
                Rural_pop_1980_capita = Rural_pop_1980/Pop_1980,
                Pwnaionwe_1980_capita = Pensioners_1980/Pop_1980,
                US = c(1,1,1,1,0),
                SB_1975_capita = SB_1975/(Pop_1975/1000),
                SB_1976_capita = SB_1976/(Pop_1975/1000),
                SB_1977_capita = SB_1977/(Pop_1975/1000),
                SB_1978_capita = SB_1978/(Pop_1975/1000),
                SB_1979_capita = SB_1979/(Pop_1975/1000),
                SB_1980_capita = SB_1980/(Pop_1980/1000),
                SB_1981_capita = SB_1981/(Pop_1980/1000),
                SB_1982_capita = SB_1982/(Pop_1980/1000),
                SB_1983_capita = SB_1983/(Pop_1980/1000),
                SB_1984_capita = SB_1984/(Pop_1980/1000))

vars <- c("Average_salary_1980",colnames(data)[c(32:46,48:ncol(data))])

means_poland <- data%>%
  dplyr::filter(., US == 0)%>%
  dplyr::summarise_at(vars, mean, na.rm = TRUE)

means_upper_silesia <- data%>%
  dplyr::filter(., US == 1)%>%
  dplyr::summarise_at(vars, mean, na.rm = TRUE)


means <- as.data.frame(cbind(t(means_upper_silesia), t(means_poland)))

means <- means[c(1,4,5,6,7,12,13,16,2,17:26),]

vars <- c("Average salary 1980 (in PLZ)",
          "Employed in industry 1975", 
          "Employed in industry 1980", 
          "Employed in agriculture 1975",
          "Employed in agriculture 1980",
          "Urban population 1975",
          "Urban population 1980", 
          "Pensioners",
          "NSZZ members 1981", 
          "SB agents 1975 (per 1,000 pop.)", 
          "SB agents 1976 (per 1,000 pop.)", 
          "SB agents 1977 (per 1,000 pop.)", 
          "SB agents 1978 (per 1,000 pop.)", 
          "SB agents 1979 (per 1,000 pop.)", 
          "SB agents 1980 (per 1,000 pop.)", 
          "SB agents 1981 (per 1,000 pop.)", 
          "SB agents 1982 (per 1,000 pop.)",
          "SB agents 1983 (per 1,000 pop.)", 
          "SB agents 1984 (per 1,000 pop.)") 

means$var <- vars
means <- means[,c(3,1,2)]
colnames(means) <- c("","Upper Silesia", "Poland")

print.xtable(xtable(means), 
             type = "latex", 
             file = "PUT YOUR FILEPATH HERE")



